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Abstract 



We present a new molecular-dynamics algorithm for integrating the equations of motion for a system 
of particles interacting with mixed continuous/impulsive forces. This method, which we call Impul- 
sive Verlet, is constructed using operator splitting techniques similar to those that have been used 
successfully to generate a variety molecular-dynamics integrators. In numerical experiments, the Im- 
pulsive Verlet method is shown to be superior to previous methods with respect to stability and energy 
conservation in long simulations. 
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I. INTRODUCTION 



Purely coUisional systems were among the first to be studied by molecular-dynamics simulationEJ. These systems 
include hard spheres or hard ellipsoids, which undergo purely elastic collisions, and square-well fluids, in which an 
attractive impulse force at a particular interparticle position is present in addition to the hard-core interactions. The 
algorithms for such systems are exact, to within roundoff error, and consist of free particle motion punctuated by 
exact resolution of the impulsive collisions — the resulting phase space trajectory is discontinuous. On the other hand, 
the vast majority of current molecular-dynamics simulations are performed on systems with continuous potentials. 
For such systems, the trajectory must be approximated using a numerical timestepping scheme such as the popular 
Verlet algorithmH. 

There exist, however, systems that are neither purely collisional, nor continuous, but are hybrids of the two. Im- 
portant-examples of such systems are the restricted primitive model (RPM) for electrolyte solutions and dipolar hard 
sphereg3. In additicm, the use of hard-core potentials with attractive continuous tails is common in pcrturbative 
treatments of liquidsQ. Since the the algorithms for simulating impulsive and continuous systems are fundamentally 
different from one another, the construction of hybrid methods for mixed systems is non-trivial and little studied. Con- 
sequently, the vast majority of studies on such systems have utilized Monte Carlo simulation techniques, eliminating 
the possibility of obtaining dynamical information. 

In this paper, we present a new method for mixed hard-core/continuous potentials, which we call the Impulsive 
Verlet algorithm. This algorithm is suitable for any continuous potential, is less likely than alternatives to miss 
collisions, and exhibits good stability and energy conservation in long time simulation. In the construction of this 
new method we have been guided by recent work iArthe use of Hamiltonian splitting methods for the development of 
efhcient and stable molecular-dynamics algorithmsda. n |— ■ 

A few ad hoc hybrid methods have been constructed for mixed (hard/soft) systemsLTtj. All of these methods are 
rather similar, in that the particles are advanced according to the continuous forces by a time step using a standard 
algorithm for continuous potentials, usually some variant of the Verlet algorithm, and the trajectories are checked for 
the existence of particle overlaps at the end or during the step. If no overlaps occur, the procedure is repeated for the 
next step. If overlaps (collisions) do occur, the system is returned to its state before the step and then is advanced 
without momentum modification by the forces to the time of collision, and the momenta are then modified according 
to the rules of elastic collision. This process is repeated until all collisions have been resolved and the end of the 
time step is reached. (One major difference between the algorithms is whether overlaps are checked only at the end 
of each step, or thcoughout the step- In the former caseS, it is possible that glancing collisions are missed during the 
dynamics.) HeyesO and Suh, et alSS apply such algorithms to the restricted primitive model for electrolytes (hard- 
sphere with embedded charges in a dielectric continuum) with some apparent success. Unfortunately, as with the 
other papers on algorithms for mixed systems, no quantitative discussion on the stability or accuracy of the algorithm 
is given, making it difficult to evaluate the quality of the methods. 

The Impulsive Verlet method is developed in the next two sections, followed by a discussion m£ certain numerical 
experiments on two model systems, comparing our scheme with the algorithm used in Suh,et alx3. 



II. SPLITTING METHODS FOR MIXED DYNAMICS 



Consider a system of N particles with instantaneous positions q — {qi, q2, Qn) in d dimensions interacting 
according to a continuous potential V^({qi}), assumed for simplicity to be spherically symmetric and pairwise additive, 
that is, 

N 

^c(q) = EE'^c(g.,), (1) 

i—l j>i 

where qij = | qj — qi |, and is any smooth function of one variable. In addition, suppose the particles to have a 
hard core of diameter cr; that is, when the distance between two particles is a an elastic collision occurs that reflects 
the momentum of each particle along the collision vector. Such a hard-sphere core can be represented formally by a 
discontinuous pair potential of the form 
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We will define the energy function of the mixed system by analogy with continuous dynamics as the sum of the 
kinetic and formal potential energies: 

H{q,p)=T{p) + VUq) + Vc{q), (3) 

where, 

N 

and 

T{p) = ^p'^M-'p, (4) 

is the kinetic energy {M is the mass matrix) and p = (j)i,p2, ...tPn), where each pi is a d-dimensional vector. Despite 
appearances, this energy function is not, properly speaking, a Hamiltonian. Nonetheless, we can view the dynamics 
of the hard-sphere fluid as the limiting dynamics in repulsive inverse-power potentials of the form Vswii") — ^/r^ , with 
/3 a large positive integer. In this sense and for the purpose of constructing numerical methods, we can interpret the 
formal energy H as representing a very hard repulsive inverse-power Hamiltonian. We will often refer to H as the 
pseudo-Hamiltonian. 

We define the flow map as the generator of the phase space trajectory, 

^^^^*^^-v...r£n. (5) 



The family of flow maps is closed under composition, 

for any times ti and ^2- n 
A continuous Hamiltonian system can often be split into integrable subproblems with Hamiltonians Hi and H2U: 

H{q,p)=Hiip) + H2{q) . (7) 

The flow map of the full Hamiltonian can then be approximated as the concatenation of flow maps for the subproblems. 
There are a variety of ways of doing this, but the most common is based on a Trotter factorization 

i'h.H = V'l.i/j ° i'h.Hi ° i'^M2 + '^(fT-^) I (8) 

where h is the time step. For a^eparable Hamiltonian such as Eq. || with Vhs — 0, this factorization reduces to the 
usual velocity-Verlet algorithmllil when Hi — T(p) and H2 = K;(q)- 

The splitting framework for continuous Hamiltonians suggests a means of constructing integrators for mixed impul- 
sive/continuous systems. A natural splitting for the pseudo-Hamiltonian is to let Hi = T(p) -f Viis(q) and H2 = Vc{q). 
(Note that in this case Hi is a function of both p and q, but since this represents a system with free particle motion 
punctuated by elastic collisions, it is exactly integrable.) This gives 

„n+l \ f Q" \ 

pn+1 j =^^,V^ O^Ph,T+V^, °^i,Va I p" ) ' 

where q" and p" are the approximations to the phase space variables after the n-th time step. In other words, the 
momenta are adjusted at the beginning of each time step by one-half step according to the continuous forces ("kick"). 
The positions are next advanced for one time step, resolving all elastic collisions, but without further momentum 
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modification by the continuous forces ("push"). At the end of the step, the momenta are advanced again by a half 
step using|-the forces calculated from the new positions (another "kick"). This is nearly identical to the algorithm of 
Suh, et a/.E3 except that there momenta are only defined at half steps and a leap-frog formulation is used: 

ri+l \ / q" \ 

^n+l/2 ^ i^h.T+VHs ° i'h.vA ^n-1/2 ■ (10) 

P /Suh VP W 

Viewing the hard-sphere potential as being approximated by a very hard inverse-power repulsive potential, we see 
that either of the above two splitting methods is symmetric (i.e. time- reversible). From Eq. |g we naively expect that 
such a method (applied to the inverse-power potential approximation) is second order accurate, meaning that in one 
step a local error of size 0{h?) is introduced; on a finite fixed time interval, these errors accumulate, but the total 
growth or global error is at most 0{h?). However, the demonstration of third-order local error requires a solution, 
and this assumption will break down in the limit of hard-sphere dynamics, in particular during a collision step. In 
fact, the local error introduced during a collision is really 0{h). 

We illustrate this point with the simple example of a nonlinear "impact oscillator" with one degree-of-freedom 
pseudo-Hamiltonian 

^^ = y + '/'hs(g) + <^c(9), (11) 

describing a point mass acted on by some potential (j)c in coUisional dynamics with a hard wall aX q = 5. The 
particle moves in the continuous potential according to Newton's equations, until an impact, when q = S, then the 
momentum changes sign and the motion continues from the impact point. 

Consider a numerical step from the point {qo,po) at time t — ioi a timestep of size h during which the particle 
motion includes a single collision event. (We mostly use subscripts to index particle number and superscripts for 
timestep, but for the discussion that follows we need to indicate powers of the momenta; so for this one-particle 
model, we will use subscripts for the timestep index.) We need to compute the local energy error contribution during 
a collisional step for this single degree-of-freedom model problem. The sequence of computations is 

?' = Po-^0;(go), (12) 

—, (13) 

K^h^h#, (14) 
p=-p, (15) 
qi = 6 ~ hp, (16) 

Pi=P-^0c(9i): (17) 

where h^f^, and h\, are the time to the next collision and the time from that collision to the end of the time step, 
respectively. 

Substituting the endpoint values into the energy relation, we quickly find 

AH = H{qi,pi) - H{qo,pa) = ^iP - ^€(^1))^ + Mn) ~ \pI " d^c[qo). 

= \{-pp + \<l>'Mp)-\d,'Mi)? 

+4>c{qi) - ^pI - 0c (go)- 

Expand (j)c in a Taylor series about q = S, substitute, and cancel like terms to obtain 

AH = -^pomqo) - €{qi)) + y (00(90) - €{qi)? 

+cj^'A5){qi - qo) + l^Umqi - sf - (go - sf) + e^. 



4 



The remainder Eh contains terms of order the third power of h or higher, i.e. \Eh/h^\ is bounded for all /i < 1 such 
that the step contains a collision. Indeed, ^{<p'ciQo) ~ is also of this order, since qi — qq is proportional to h. 

From the equations 

qi = 5 - hbP, qa = S ~ h^p, 

and the use of a Taylor series expansion of 0^, we arrive after discarding terms of order three or higher at, 

AH = hih*^cj^'^iS)p' - 4>'^[5){h, - h#)p + i^i' - hl)p^ + Eh. 

with Eh again of third order. This finally leads to 

AH = ih*-^(h - {h# + K))4>'l{5)f - 4>'A5){h - h#)p + Eh 
^ ~(t>'^{S){h ~ h#)p + Eh. 
Therefore, the expected energy error introduced in this one collisional step is 

H{q,,pi) - H{qo,po) = -(t>'ciS){h - h#)p + Eh, (18) 

where the quantity \Eh/h^\ is bounded independent of h. (A similar result would hold for the solution error.) 

Technically speaking, it is incorrect to say that the energy jump in one step is 0{h) since if we decrease the 
timestep h sufficiently, there will be no collision event within the particular step, and so the error will revert to 
0{h^). Nonetheless, in any timestepping simulation in which there are collision events, these steps will introduce 
errors proportional to h. If we define the local approximation error eioc as the maximum of magnitudes of the local 
errors introduced, then eioc is of first order in h, not third order as we would expect in the continuous case. Since there 
are, in general, a finite number of such collisions in any finite interval, the accumulation is bounded and the global 
error is also 0{h). The apparent contradiction of an odd-order symmetric method is just one of several anomalies 
that result from the complex transition from the smooth problem to the discontinuous limit. In another terminology, 
we could say that the splitting method undergoes an order reduction for stiff potential wells. 

From this discussion and Eq. we expect the naive splitting method to give rather poor energy conservation, 
except in three special cases: 

Case 1 Collisions do not occur within timesteps but precisely at the timesteps, so third order is recovered. 

Case 2 The collisions occur at precisely the middle of a timestep, so that the first order term in the error formula 
vanishes and third order local energy drift is again recovered. 

Case 3 Third order will be recovered if the derivative of the continuous pair potential vanishes for two spheres in 
contact. 

To illustrate this last point, we apply the method to one degree-of-freedom anharmonic "impact oscillator" with a 
continuous potential, 4>c{q) = \q^ + \q^. We show in Fig. |l|the maximum total energy error as a function of the time 
step when the wall is placed at g = 0.00 and q = —4.00. One can see that the naive splitting is a second order method 
when the derivative at the wall vanishes . 

Because it is only applicable for a relatively limited class of potentials the naive splitting method is not a candidate 
for a viable general technique, however, it does provide a good starting point for the development of a general method, 
which we call the Impulsive Verlet (IV) algorithm. 

III. IMPULSIVE VERLET 

To develop our method, we deliberately exploit two of the special cases in the naive algorithm for which third order 
can be expected, namely Cases 1 and 3 mentioned at the end of the previous section. (Case 2, the situation that 
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collisions occur at the midpoint of the time interval, does not appear to be of practical use.) We begin by introducing 
an artificial splitting of the continuous potential, 4>c{'lij), into into a short-ranged part, (j3i{qij), and a long-range part, 
02 according to 

0c(9ij) = 01 (gy) + 02(9-y) • (19) 

(This decomposition is similar to that invoked in multiple timesteppingBEl^Q molecular-dynamics algorithms.) For 
the reasons discussed above, the long-range (and therefore most expensive to calculate) part of the potential is defined 
so that the derivative vanishes at the hard-core separation. We define 02 (9y ) as follows: 

^■(91), q < qi, 

Mq) = { Piq): qi<q< 92, (20) 
0c(9), g > 92 , 

where, qi and 52 are parameters, and P{r) = Ao + Air + A2r^ + A^r^ is a Hermite interpolant introduced so that the 
two potentials are smooth to the order for any continuous potential. From Eqns. |l9|) and [20| ), 4>i{q) is given by 

Mq)-P{qi) q<qi 
Mnj) = { Mq)-Piq) qi<q<q2, (21) 
q>q2- 

The continuity condition, P{r2) = 4>c{'r2) , and the smoothness conditions, P'{q2) — 'P'ciq^): P'{qi) — Oj and P"(gi) — 
0, allow us to calculate the coefficients of the Hermite interpolant, giving 

^3 = ^TzV^ 2T for gi 7^ 92 , (22) 

Qri{qi-q2)+i{qi-q() 

A2 = -3giA3, (23) 

Ai = Sg^As, (24) 

Ao = -{Aiq2 + A2ql + A^ql) + Vc{q2) ■ (25) 

(An example of this potential splitting for an inverse-sixth-power attractive potential, 0c('?) — ~^(s'/q)^, is shown in 
Fig.|.) 

Next, we define A^-body potentials Vi and V2 as a sum of pair contributions from 0i and 02, respectively. We then 
split the total Hamiltonian in the following way: 

Hi{q,p)^T{p) + VUq) + Vi{q) (26) 

and 

H2{q) = V2{q) . (27) 

The Trotter factorization (Eq. ^) is now applied to this splitting. The problem now is that H2 is not integrable and 
its flow map must be approximated. This is done is the following way: 

"c + l 

i=l 

where ric is the number of hard-sphere collisions between during the time step h, t^'^^ is the time between each collision 
(with T-[^'' being measured from the beginning of the time step until the first collision and t^^^ measured from the 
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last collision to the end of the time step), and ipVi^a is an operator representing the resolution of each elastic collisions. 
This is essentially the execution of a Verlet step of length t^'^'> between each elastic collision. The collision times can 
be calculated since the Verlet step generates a quadratic trajectory, which together with the collision condition for 
two particles i and j can be written as 

||g,(r(^))-9,.(r(^))f -a2=0, (29) 

generates a quartic equation for t^'^\ 

We describe below the algorithm for the Impulsive Verlet molecular-dynamics simulation in more detail. 



Impulsive Verlet Timestepping Algorithm 
pr''"'°=P"'° + iF2,,(q-0)/l 

do ic = 1, Tic 

pn+l/2,>-l/2 ^ p„+l/2,,-l ^ lp^_^(qn,»-l)^i 
qn,i _ q7i,i — 1 _j_ ■|y£-lpn + l/2,v^j 
^n+l/2a ^ p„+l/2,>-l/2 _^ lFj^^(qn-)^i 

_,n + l/2,i _ ,/, { q"'*'' A 
P - VVhS y ^n + l/2,i J 

end do 

pn + l/2,n, + l/2 ^p„ + l/2,„, ^ i Fi (q"."^ " Er=l ) 
qn + 1,0 ^ qn,n, _^ lpn + l/2,ne + l/2 _ ^"c^ 

pn + l/2,n, + l ^ p„ + l/2,„, + l/2 _^ 1 F^ (q" + l .0 ) (/l - J2Zl 
pn + 1.0 ^ pn + l/2.„<= + iF2,,(q" + 1.0)h 



To make sure that no collisions are missed it is necessary to ensure that the quartic equation (Eq. |2^) is accurately 
solved to give the nearest root to zero. This is not a trivial problem as the solution becomes increasingly unstable as 
smaller time steps are used (i.e. when the time to collision is small) . To ensure the inaccuracies are not large enough 
to affect the overall accuracy and order of the method, we employ Laguerre's methodtj to find all roots of the quartic 
and take the smallest, positive real root, which is then refined using Newton-Raphson. This proved to be sufficient 
at all but the very smallest time steps studied. 

There is a small probability that the Impulsive Verlet method can miss a grazing collision, since the trajectories 
that are followed in determining collisions are quadratic approximations. However, this probability is greatly reduced 
in comparison to the method of Suh, et al. or any other algorithm that uses linear motion to determine the collisions. 

IV. NUMERICAL EXPERIMENTS 

We test the Impulsive Verlet algorithm using as our continuous potentials, 4>c{q), the Lennard- Jones potential 



4e 



12 / \ 6 

a \ I a 



(30) 



and an attractive inverse-sixth-power potential 



,9 



(31) 
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In both potentials a is the same as the hard-core diameter. We truncate both potentials at the distance q* = qc/cr — 2.5 
and, to ensure their continuity, they are shifted so that the value of the potential at the cutoff is zero. In implementing 
the Impulsive Verlet algorithm, we split each potential as prescribed in Eq. |l^ |2^, with qi and (72 as input parameters. 
For theXennard- Jones potential there is, of course, a natural splitting, namely that of Weeks, Chander and Anderson 
(WCA)a, where the potential is split at the minimum with ql = q2 = 2^/^, which gives the following splitting: 

</)i,„(g;WCA) = ^'^^^ ^ ^ ^ ' (32) 

I 0, 5 > 26(7. 

'^2.Lj('?, wcA) = 1 4^[(^)_i;': (^)-6]^ 1 1 X\a, (33) 

The MD simulations were carried out on systems of 108 particles. The system of reduced units was chosen so that all 
quantities are dimcnsionless. So, as units of distance and energy we used the potential parameters g and e, respectively, 
and the mass of one atom as the unit mass. The unit of time is (mtr^/e)^^^. An asterisk superscript indicates reduced 
units. Except were otherwise indicated all simulations are performed using a reduced density p* — pa^ = 0.9 and 
reduced temperature T* = kT/e — 2.5. In addition, a cubic box jSsdth periodic boundary conditions is used. For 
greater efficiency, the MD program incorporates three neighbor listsES for the evaluation of the short-range force, the 
long-range force, and the collision times. 

The results of the Impulsive Verlet on the instantaneous total energy for the Lennard-Jones and the attractive 
inverse sixiii continuous potentials are illustrated in Fig. ^ and |^. A comparison to the naive splitting algorithm of 
Suh, et alS3 is also made for both potentials. The superiority in energy conservation and stability of the Impulsive 
Verlet algorithm over the naive splitting method is striking. 

We study in Fig. ^ the order of the method while varying qi and 52 ■ The order is obtained by plotting (on a log- log) 
the maximum energy error for a fixed-length simulation versus the time step. A comparison with a straight line of 
slope two tells us that the method is of second order for various values of ql and q2 ■ (Note the slight deviation of the 
slope at very small time steps from the theoretical value of 2.0 is due to the difficulty in solving the quartic equation 
for the collision times when the time to collision is very small. This is not a real problem in practice since the goal of 
molecular-dynamics simulation is to use the largest time steps possible.) 

Finally, to demonstrate the ability of the Impulsive- Verlet method to yield relevant dynamical quantities, we show 
in Fig. ^ the result for the normalized velocity autocorrelation function , C{t) = (v(t) • v(0))/(v(0) • v(0)), for the 
Lennard-Jones system (108 particles) with p* = 0.9 and T* = 0.9. In this calculation we use a splitting with ql — 1.122 
and ^2 = 1-5- 

V. CONCLUSION 

We have introduced a molecular-dynamics method for mixed hard-core/continuous potentials, which we refer to as 
the Impulsive Verlet algorithm. This algorithm is produced by extending general potential splitting methods to the 
specific case of mixed potentials. In addition to providing a mechanism for generating the Impulsive Verlet method, 
the potential splitting formalism helps to understand the failings of previous methods. The Impulsive Verlet algorithm 
uses a quadratic trajectory between collisions and does not miss any collisions of the approximate trajectory. As a 
result the algorithm is suitable for any type of continuous potential, is second order, has good energy preservation, 
and is far more stable over long time simulation than previously integrators for such systems. (A detailed theoretical 
analysis of the algorithm is the subject of current research.) 
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FIGURES 



FIG. 1. The maximum total energy error as a function of the time step for one degree-of- freedom anharmonic 

"impact oscillator" is using the naive splitting approach. The wall is placed at q = —4.00 (square) and at q = 0.00 
(circle). Comparison is made with lines of slope two(solid line) and one(dashed line). 

FIG. 2. A potential splitting of a inverse-sixth-power attractive interaction, — (^)^, with qi/a = 1.1 and 
92 /c = 1.200. The short range and the long potentials are (a) and (b) respectively. 

FIG. 3. Instantaneous total energy for a 108 particle simulation using a Lennard-Jones continuous potential with 
a hard-sphere core, using both Impulsive Verlet (solid line) and the naive splitting algorithm of Suh, et al. (dashed 
line). The time step is /i* = 4 x 10~^. 



FIG. 4. Instantaneous total energy for a 108 particle system interacting via an inverse sixth-power attractive 
potential with a hard-sphere core, using both Impulsive Verlet(solid line) and the naive splitting algorithm of Suh, et 
al. (dashed hue). The time step is /i* = 4 x 10""^ . 



FIG. 5. The maximum total energy error as a function of the time step for a system of 108 particles using the 
Impulsive Verlet algorithm and the Lennard-Jones potential. Comparison is made with a line of slope two. 



FIG. 6. Normalized velocity autocorrelation as a function of time for 108 Lennard-Jones particles at p* = 0.9 and 
T* = 0.9, using the Impulsive Verlet algorithm with a time step h* = 1 x 10~^. 
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